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ABSTRACT 

The helicity condensation model has been proposed by Antiochos (2013) to explain the 
observed smoothness of coronal loops and the observed buildup of magnetic shear at filament 
channels. The basic hypothesis of the model is that magnetic reconnection in the corona causes 
the magnetic stress injected by photospheric motions to collect only at those special locations 
where prominences form. In this work we present the first detailed quantitative MHD 
simulations of the recon n ection evolution proposed by the helicity condensation model. We use 
the well-known ansatz of modeling the closed corona as an initially uniform field between two 
horizontal photospheric plates. The system is driven by applying photospheric rotational flows 
that inject magnetic helicity into the system. The flows are confined to a finite region on the 
photosphere so as to mimic the finite flux system of, for example, a bipolar active region. The 
calculations demonstrate that, contrary to common belief, coronal loops having opposite helicity 
do not reconnect, whereas loops having the same sense of helicity do reconnect. Furthermore, we 
find that for a given amount of helicity injected into the corona, the evolution of the magnetic 
shear is insensitive to whether the pattern of driving photospheric motions is fixed or quasi- 
random. In all cases, the shear propagates via reconnection to the boundary of the flow region 
while the total magnetic helicity is conserved, as predicted by the model. We discuss the 
implications of our results for solar observations and for future, more realistic simulations of the 
helicity condensation process. 
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1. INTRODUCTION 


High spatial resolution observations of the Sun have revealed at least two major, puzzling 
features of the eorona. XUV and X-ray images have shown that eoronal loops are smooth and 
laminar, with little or no tangling (Schrijver et al. 1997). This result is difficult to understand, 
given that the quasi-random motions of the photosphere continuously drive the field. In fact, the 
tangling of the field by photospheric motions is widely believed to be essential for the observed 
coronal heating (Parker 1988; Klimchuk 2006). Another long-standing puzzle is the magnetic 
shear in filament cha nn els (Martin et al. 1992; Zirker et al. 1997; Martin 1998; Mackay & van 
Ballegooijen 2001; Pevtsov et al. 2003a). Multi-wavelength observations ranging from white 
light to X-ray have shown that wherever the normal flux at the photosphere changes sign, at so- 
called polarity inversion lines (PILs), strongly sheared field inevitably builds up in the corona 
leading to the well-known phenomenon of prominences and filaments (Priest 1988; Tandberg- 
Hanssen 1995). 

In recent work, Antiochos (2013) proposed that a single process, helicity condensation, 
explains both the laminar coronal loops and the filament channel shear. Magnetic helicity is the 
topological measure of the total linkages of the field lines in any flux system, and has been 
defined rigorously for a coronal domain (Berger & Field 1984; Finn & Antonsen 1985). The 
arguments leading to the helicity condensation model are straightforward. First, the small-scale 
granular and supergranular flows at the photosphere must continually stress the coronal magnetic 
field, building up small-scale current sheets that inevitably dissipate by reconnection. This is the 
basic scenario behind the so-called nanoflare model for coronal heating (Parker 1988; Klimchuk 
2006): the free energy injected into the corona by photospheric motions is continuously 
dissipated by impulsive reconnection events. Second, if these photospheric motions also inject a 
net helicity into the coronal field, then the helicity can only accumulate because it is conserved 
under reco nn ection (Woltjer 1958; Taylor 1974, 1986; Berger 1984). In fact, a large number of 
observations ranging from subsurface flows (Duvall & Gizon 2000; Gizon & Duvall 2003; 
Komm et al. 2007), to sunspot whorls (Pevtsov et al. 2003b), to filament channel chirality 
(Martin et al. 1992; Zirker et al. 1997; Martin 1998; Mackay & van Ballegooijen 2001; Pevtsov 
et al. 2003a) show a predominant net helicity in each solar hemisphere, with negative-helicity 
(left-handed) structures in the north and positive- (right-handed) in the south. The puzzle. 
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however, is that coronal loops do not show the tangling or twisting that is expected due to this 
helicity buildup. 

Third, in systems evolving via turbulent reconnection, the helicity undergoes an inverse 
cascade in phase space and “condenses” at the largest system wavelength (Biskamp 1993). In the 
case of the corona, the relevant largest wavelength is the scale of the whole flux system, which is 
defined by the PIL (Antiochos 2013). As a result of this helicity condensation process, the 
magnetic flux tubes far from the PIL remain relatively unstructured, whereas the flux near the 
PIL becomes strongly sheared. This is exactly what is needed to account for both laminar loops 
and filament channels. Furthermore, Mackay et al. (2014) incorporated helicity condensation into 
large-scale magneto-frictional simulations of filament-channel formation, and showed that the 
process can reconcile major discrepancies between observations and previous theories in the 
dominant skew of polar-crown filaments. 

In this paper, we investigate several critical assumptions underlying the helicity condensation 
model using full MHD simulations, described in Section 2, that capture rigorously the injection 
and reconnection dynamics. First, we test in Section 3 the conjecture presented by Antiochos 
(2013) that only flux tubes with the same sign of helicity reconnect, and that tubes with opposite 
helicity do not. This result seems highly counterintuitive, but if correct, it places important 
constraints on coronal evolution: it implies that helicity cannot simply be canceled by 
reconnection. Second, we investigate in Section 4 the effect of the pattern of photospheric flows 
on helicity condensation. The results derived in Antiochos (2013) were based on modeling the 
photospheric flows as spatially fixed rotations, but observed granular and supergranular flows 
are clearly much more complex, appearing and disappearing quasi-randomly on the photosphere 
(Rieutord & Rincon 2010). The helicity condensation model is robust if the process is insensitive 
to whether the pattern of flows is stationary or time-varying; the only requirement is that the 
flows must inject a consistent handedness of helicity over an extended region large compared to 
the individual rotations. To test the model, we simulate helicity injection and evolution using 
both fixed and random patterns of photospheric flows. The results are compared to determine, in 
particular, whether the helicity condenses along the PIL as predicted, irrespective of temporal 
variations in the flow pattern. Finally, we discuss the implications of our results for solar 
observations and for future theoretical and modeling studies. 
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Figure 1. Schematic of the simulation system consisting of a uniform vertical magnetic field 
between two horizontal, photospheric plates. The system is stressed by applying rotational 
motions at the photosphere, indicated by the circles and arrows. Left panel: co-helicity case; 
right panel: counter-helicity case. 


2. NUMERICAL MODEL 


In this exploratory numerical investigation of magnetic helicity injection and evolution, we 
adopt the simplest possible coronal model consisting of an initially uniform magnetic field with a 
gravity-free Cartesian atmosphere bounded by planar “photospheres” at both the top and bottom 
of the domain. For such a system, a coronal “loop” corresponds to a straight flux tube extending 
from the bottom plane to the top plane of the domain. The initial field is vertical and the field 
strength, plasma density, and temperature are uniform throughout. This model was originally 
proposed by Parker (1972) for investigating the effects of photospheric motions on coronal 
heating and has been used by a number of authors in numerical studies of the corona (Mikic et al. 
1989). At both bounding photospheric planes, we impose horizontal cyclonic flows that are 
incompressible and so leave the normal component of the magnetic field undisturbed. These 
simple flows minimize the complexity of the simulated evolution, but capture the properties of 
the solar surface convection that are critical to the injection of magnetic helicity into the corona. 

Figure 1 is a schematic diagram of the simulation setup that we use. The twisting flows are 
imposed at both planes simultaneously, with the sense of rotation reversed from one plane to the 
other as viewed from outside of the Cartesian domain, to maximize the helicity injection. For 
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clarity, we show only two rotations, but given sufficient computational resources one could 
impose an arbitrary number. From the perspective of the figure, we are looking down upon the 
“photosphere” that is the bottom plane, but up from below through the “photosphere” that is the 
top plane. Thus, as seen from within the corona - e.g., from the center of the domain, at the 
“apex” of our “loop” - the sense of the cyclonic motions is the same at both footpoints, inducing 
maximum twist within the loop. 

We calculate the resultant evolution of the coronal magnetic field and plasma by solving the 
magnetohydrodynamics (MHD) equations in three Cartesian dimensions and time in the form 

^ + V.(pv) = 0, (1) 

^+V»(pvv) = -VP+T(VxB)xB, (2) 

dT 

— +V. rv)=(2-r)rv.v, (3) 

ot 

— -Vx(vxB) = 0. (4) 

All symbols have their usual meanings: /?is mass density, T is temperature, P is thermal 
pressure, y is the ratio of specific heats, v is velocity, B is magnetic field, and t is time. The 
ideal-gas equation of state relates the first three variables through the gas constant R, 

P = pRT. (5) 

Our focus in this paper is on understanding the magnetic-helicity injection, accumulation, and 
transport, not the details of the thermal response of the coronal plasma, so we adopt the adiabatic 
description of the plasma energy evolution. 

We solve the prescribed time-dependent equations of MHD with the Adaptively Refined 
Magnetohydrodynamics Solver, ARMS (DeVore & Antiochos 2008), which employs Flux- 
Corrected Transport algorithms (DeVore 1991) and a finite-volume representation of the 
variables to obtain its solutions. The integration methods introduce necessary, stabilizing 
numerical dissipation terms into the nominally ideal equations written above. These terms are 
negligibly small in regions where the solution is smooth and well resolved on the grid, but are 
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large enough to keep the solution stable and monotone in regions where any variables change 
substantially on the grid scale. At electric current sheets associated with discontinuities in the 
direction of the magnetic field, the numerical dissipation term in the induction equation allows 
reconnection to occur. This will be evident in the results presented below. 

Due to the homogeneity of the MHD equations, scale factors for mass, length, and time - or 
any equivalent combination of those three fundamental scales - can be extracted from the 
variables and removed from the equations. This means that any solution to the equations that is 
obtained actually represents an entire family of solutions that adhere to the underlying scaling 
relations. We exploit this freedom to extract scale factors Us, Ps, and Bg for length, mass density, 
and magnetic field strength, respectively, from the equations. The ideal gas law allows us also to 
extract a temperature scale factor Tg. 

After performing this rescaling, we prescribe the initial state of our system to be 

p(x,0) = Po, 

t(x,o) = t;, 

P(X,0) = P„ (6) 
v(x,0) = 0 , 

B(x, 0) = BqX. 

The vertical direction of our domain lies along the x coordinate. We set the constants in the 
initial state, together with the radius ao of our cyclonic convection cells (see below), to the values 

= .125, 

A =1’ 

To=h (7) 

Po=.05, 

B^ = 

These choices fix the values of the gas constant, R = 0.05; the initial uniform Alfven speed, cao = 
Bo^ATipf) = 1 ; and the initial uniform plasma beta (ratio of thermal to magnetic pressure), Pq = 
SnPo/Bo =0.1. The last selection ensures that the evolution is magnetically dominated, as is 
generally true in the corona. In addition, our simulation domains all have unit height, so AlfVen 
waves in the initial configuration require one unit of time to propagate between the top and 
bottom of the domain; i.e., the simulation time is expressed in Alfven transit times. 
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We impose open-boundary, zero-gradient conditions at all four side walls (the y and z 

bounding planes) of our domain, on all variables. That is, 

dp dr dv dB 

— = — = 0 , — = — = 0 , ( 8 ) 
dn dn dn dn 

where n =y,z is the normal coordinate. We apply the same conditions at the top and bottom (the 
X bounding planes) of our domain, with the exception of the velocity, where we prescribe 

V, =xxV;if(y,z,t) (9) 

for the normal (v„) and tangential (vt) components, respectively. This incompressible flow 
pattern leaves the values of the scalars (p and T) and the normal magnetic field (Bn) undisturbed 
at our photospheric planes, but induces a tangential component of the magnetic field (Bt) 
wherever the stream function % has a finite gradient. We impose a superposition of simple, 
separated, cylindrically symmetric flows centered at positions (yo,zo) and having finite radial 
extent r = (y~yo) +(z-zq) < ao , taking the form 

x{y,z,t) = VnaJ(t)g(r), (10) 


where y(t) is a prescribed temporal profile (given below) and the spatial profile g(r) is 


g{r) = ■ 
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Beyond r = ao, we fix g(r) = 0. The corresponding angular rotation rate is 


Q(r,t) = v^ajit)-^ = ~^f(t) 
r dr Uq 
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( 12 ) 


which vanishes both at the center of the circle (r = 0) and at the edge of the circle (r = ao). For 
large / and m, the fall-off near both the center and edge is steep, so that the maximum linear 
velocity |v|max in the annular region of finite flow approaches the profile constant |vo|. As stated 
above, we set ao = .125, and we also assumed moderate values / = m = 2 for the profile indices. 
The maximum linear and angular speeds of the flow patterns are found to be 

Hmax ~ l^o|2^(2/ + l) j (^2m + 2l + 1^ = .213|vq|, (13) 
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We chose vo = +0.938, so that |v|max = 0.200; this value is sufficiently sub-Alfvenic to allow the 
field to evolve approximately quasi-statically, but large enough to yield amply twisted magnetic 
fields in modest-duration simulation runs. The corresponding maximum rate of the clockwise 
rotation, |Q|max = 1.876, corresponds to a full rotation period 27i/|Q|max = 3.35. Finally, to ramp 
the flows smoothly up from zero and/or back down to zero, we used a cosine temporal profile 


/(o=i 


1-cos 


V j 


(15) 


over intervals t e \t\,t 2 \. Outside of such intervals, we fix/(0 = 0 or 1, as desired. For example, 
to obtain one half-turn of twist at each footpoint over a full ramp-up/ramp-down cycle beginning 
at some time t\, we simply set to = h and = h + 27i/|Q|max = h + 3.35. 

We used the adaptive mesh refinement capabilities of ARMS to resolve very finely the 
volume of the domain where the photospheric flows were imposed and where the coronal flux 
tubes were twisted, with allowance toward the sides to accommodate flux tube expansion. The 
mesh was composed of 8x8x8 blocks of uniform, cubic grid cells. We used 4 blocks, or 32 grid 
points, to span the radius «o of each of our cyclonic convection cells. The full convection 
pattern of multiple cyclonic cells, the narrow lanes between them, and a buffer zone around the 
perimeter of the pattern were covered by fine grid blocks. This finely gridded region extended 
vertically throughout the corona. Outside of this region, the grid was allowed to coarsen by two 
levels toward the side walls, with the spacing increasing by a factor of two at each refinement 
change. The resulting grids were about 50% smaller than uniformly refined grids. 

To monitor and diagnose the simulations, we use five global quantities: the total kinetic {K) 
and magnetic {M) energies and magnetic helicity {H) in the volume, and the average magnetic 
flux densities through the vertical y = 0 {Fy) and z = 0 (Fz) planes. The energies and flux 
densities were computed directly as volume and area integrals over the grid during the runs, 

K=Apv^dV, ( 16 ) 

V ^ 

M=\—B^dV, (17) 
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(18) 



We calculated the magnetic helicity both as a volume integral of the instantaneous value (Finn & 
Antonsen 1985) and as a surface integral of the corresponding injection rate, 

// = j(A + A,).(B-B,)rfF, (20) 

V 

^ = + (21) 


In these formulae, B is the instantaneous magnetic field in the system with vector potential A, 
and Bjp is the potential field with the same normal component at the boundary as B. Note that in 
our case Bp is constant in time and is simply the initial uniform field. Both integrals were 
evaluated numerically in a post-processing step using the ARMS simulation data; the surface 
integral also can be performed analytically. We chose for the vector potential Ap of our uniform 
initial potential field the symmetric form 

( 22 ) 

then Bp = V xAp = 5 qX , and Ap«Bp = 0. The instantaneous vector potential at later times was 
computed by integration on the adapted grid, using Ap as the boundary value at our bottom 
photosphere at x = 0: 

X X 

A (x, y , z, t (y, z) -h y J dxB^ (x', y, z, t ) - z J dx'B^ (x', y,z,t). (23) 

0 0 


This value is then used in the volume integral for H. The surface integral for the helicity 
injection rate dHJdt over a cyclonic convection cell centered at (yo,zo) is evaluated by noting that 


A(y,z) = A(yo,Zo) + ^[(y-yo)z-(z-zJy] 
= A(yo,Zo) + Y^9. 


(24) 
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where cp is the cylindrical angular coordinate about the cell center, whence 

^ = -2Sj(A.v)rfS 

= -25q A (jo , Zq ) • j vJ5' - Bl j r(p • \dS 

Sc s, 

= -Bl\r'^Q.{r)dS (25) 

•5c 

«0 

= -ItiBI j r^Cl{r)dr 
0 

^ ^0 7^2 4 

— +71 — — , 

(/ + 2)(m + / + 2) 


where we have used the fact that our prescribed vector velocity v averages to zero over each cell 
and substituted its explicit radial dependence. The average helicity injection rate is one-half of 
the above peak rate over a ramping cycle for the flows of duration r, during which the helicity 
injected is 


AH„ = = +7T 


m 


dt 2 [l + 2){m + l + 2) 2 


^0 ^ 

OqUq . 


(26) 


For our parameters, these expressions yield dHJdt = 6.0x10 ^ for the peak rate, and A/fc = 
1.0x10 for the total helicity injected over a full cycle of duration r= 27i:/|Q|max = 3.35, by each 
rotation comprising the flow pattern. We note for completeness that at our top photosphere the 
vector potential in general will differ from its initial value used in the helicity surface integral 
above, but only by the horizontal gradient of some scalar function y/{y,z), since the normal 
magnetic field B„ is unchanged. The helicity injection rate also is unchanged, as we intuitively 
expect, because the additive term vanishes: 

J(^V^//'*v)t/5' = Jv* (^//'v) = I y/\ • rdL = 0, (27) 

Sc s^ 4 


since v is incompressible and tangential to (indeed, vanishes at) its perimeter curve Lc. Finally, 
the sign of B*h reverses at the top boundary (n is the normal direction into the coronal 
volume), but we also reverse the sign of vo (and thus v) there to compensate; consequently, the 
helicity injection rate has the same sign and magnitude at the top and bottom photospheres for 
each of our cyclonic convection cells. 
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3. INTERACTIONS OF TWO TWISTED FLUX TUBES 


To gain physical insight into the helicity condensation process, and to test the conjecture in 
Antiochos (2013) that only like-helicity loops reconnect, we first investigate the most basic 
configuration having only two rotations at each of our photospheres (Figure 1). This flow pattern 
produces two twisted flux tubes in the corona of our numerical model. The domain is the unit 
cube [0,l]x[-0.5,+0.5]x[-0.5,+0.5] which, if uniformly gridded, would have been spanned by a 
256^ grid. Paired flow circles at the bottom (x = 0) and top (x = 1) planes are positioned at 
identical horizontal centers (yc, Zc), which for these two-flow cases are set at yc = ±ao, Zc = 0, 
where ao = 0.125. Thus, the vertical y = 0 plane separates our initially straight flux tubes that 
become twisted, with the outer edges of the two tubes just touching at (y,z) = (0,0) for all xatt = 
0. Note that the plane y = 0 is a symmetry plane for the given domain and imposed rotational 
motions, irrespective of the sense of the motions (see Figure 1). This condition is very useful for 
analyzing the resulting evolution, in particular for determining exactly when and how much 
reconnection occurs. 

In these simulations, the twisting flows were ramped up to full speed over time t e [0.0,0. 5], 
using to = 0 = 0 and t 2 = 1 in (15), and then held fixed thereafter. For small twist, i.e., at early 
times t in the simulation, the imposed footpoint motions create a single flux tube in each of the 
two half-volumes, y < 0 and y > 0. Therefore, we expect to find 5^; = 0 at y = 0 to a very high 
degree of approximation, with no magnetic flux crossing the plane y = 0 or linking the flux tubes 
to each other. On the other hand, for sufficiently large twist, at late times t, the two tubes will 
expand to come into contact and interact. If they reconnect, then linkages will form across the y 
= 0 plane; in other words, the twist flux injected into each axial flux tube separately will merge, 
so that it encompasses both tubes. If there is no reconnection, it will eventually become 
energetically favorable for the tubes to kink by bending around each other and becoming helical. 
Such an ideal evolution also will break the symmetry about the y = 0 plane, but will not establish 
any linkages between the tubes. In either case, the evolution is expected to develop nonzero By at 
y = 0 at some point in the simulation; that this indeed occurs is shown below. 
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3.1 Co-Helicity Case 


Let us first consider the interaction of two flux tubes with the same sense of twist, i.e., sign of 
helicity. We apply the boundary flows shown in the left panel of Figure 1, counterclockwise at 
the top plane (x = 1) and clockwise at the bottom (x = 0). Because the flows are identical for the 
two flux tubes, the twist component in each is identical, at least during the early ideal evolution. 
As is evident from the left panel of Figure 1, this implies that the twist components are 
oppositely directed at the contact line between the two flux tubes, (y,z) = (0,0), so that a large 
current builds up there. Alternatively, we note that stagnation points occur in the flow at their 
contact points, y = 0, on the boundary planes, which implies that an exponentially growing 
current will form in the corona along the (y,z) = (0,0) line (Antiochos & Dahlburg 1997). The 
growth of the currents eventually will lead to reconnection, but only the twist components of the 
field in the flux tubes can reconnect; the axial components are always parallel, as is evident from 
Figure 1. 

Figure 2 shows the evolution of six selected field lines in response to the boundary driving. 
Results are shown at five representative times during the simulation {t = 3, 4.25, 4.5, 6, and 10). 
The field lines in each frame of the figure are traced from the footpoints on the x = 0 plane along 
the z = 0 line, at positions y = ±0.005, ±0.0625, and ±0.245. These positions also determine the 
colors of the lines in each frame and are marked by the colored beads; note, however, that the 
field lines are traced from fixed grid positions, not from fixed footpoints that convect with the 
flow. Consequently, the field lines in the various frames are not the same, even if they have the 
same color. 

The left column shows magnetic field lines and contours of plasma velocity magnitude on the 
X = 0 plane viewed from above. The right column shows a side view of the two flux tubes in 
which the X = 0 (bottom) and x = 0.9 planes are drawn with contours of velocity magnitude on 
both planes. At t = 3 and t = 4.25, the magnetic field lines in each flux tube are clearly twisted, 
but confined to their individual flux tubes. Near the y = 0 plane, the field lines (white and black) 
from the y > 0 and y < 0 domains have already developed opposite tangential components, 
however there is not yet any evidence for reconnection. 
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Figure 2. Evolution of six field lines at seleeted times t during the co-helicity simulation. Left: x 
= 0 plane viewed from above; right: x = 0 (bottom) and x = 0.9 planes. Velocity magnitude is 
color-contoured on the planes. 
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By time t = 4.5 (At = 0.25 later), the dark-blue field line rooted in the left-side flux tube has 
crossed the y = 0 plane to circle the right flux tube and, symmetrically, the light-blue field line in 
the right-side flux tube has crossed that plane to surround the left tube. The penetration of these 
field lines through the y = 0 plane is a definitive indicator of magnetic reconnection. This 
reconnection is exactly the type of evolution required by the helicity condensation model 
(Antiochos 2013). By time t = 6, all the field lines at the selected positions have reconnected and 
crossed from one flux tube to the other. Finally, by time t= 10, we observe a clear trend for the 
field lines to twist around the outer boundaries of the two flux tubes. Note especially the dark- 
and light-blue lines, which nearly encircle the two tubes. 

An effective procedure for analyzing the reconnection evolution of the two flux tubes is to 
examine the behavior of the magnetic field component By. As described above, any linkage 
between the two flux tubes that forms as a result of reconnection will be seen immediately as the 
appearance of By flux crossing they = 0 plane. Figure 3 displays the By and B^ components on the 
horizontal plane one quarter of the way above the bottom boundary, at x = 0.25, at the simulation 
times t = 3, 4.25, 5.5, 6, and 10. At the early times t = 3 and 4.25, we see the signatures expected 
for two separate twisted flux tubes. There is no apparent By flux crossing they = 0 plane and B^ is 
well-organized and symmetrical. Furthermore, at t = 3, the By contours in the left panel of Figure 
3 can be obtained simply by rotating the B^ contours in the right panel clockwise by 90°, as 
would be expected for an axisymmetric twisted flux tube. At t = 5.5, well after reconnection has 
set in (at t = 4.5, as shown in Figure 2), a small portion of the By flux contour now crosses they = 
0 plane, and the B^ flux starts becoming less organized. By time t = 6, the two same-signed By 
patterns that earlier belonged to the two separate flux tubes have coalesced, forming a complete 
circle with a single positive By pattern on one side and negative on the other, indicating that the 
two flux tubes have merged. Meanwhile, the B^ distribution also shows a dramatic change in 
Figure 3. At the outer boundaries of the two flux tubes along the y direction, B^ increases in 
strength; at the interface region between the two tubes, near y = 0, the B^ patterns fade and then 
disappear by time t = 10. The final pattern of B^ again resembles that of By rotated by 90°, but 
now for a single twisted flux tube rather than two tubes. 

The evolution of the magnetic field is even clearer from side views. Figure 4 shows By 
through the y = 0 plane (left column) and B^ through the z = 0 plane (right column), at t = 3, 4.25, 
5.5, 6, and 10. We note that at t = 3, there is no flux through they = 0 plane, and the flux through 
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Figure 3. Left: color contours of By on the x = 0.25 horizontal plane at selected times t during 


the co-helicity simulation. Right: contours of on the same plane. 
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t=3 



Figure 4. Left: color contours of By on the y = 0 vertical plane at selected times t during the co- 
helicity simulation. Right: contours of on the z = 0 vertical plane. 
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the z = 0 plane appears relatively straight and unstructured. This is exactly what is expected for a 
force-free vertical twisted flux tube, which should have Bz = constant along the axial direction. 
By t = 5.5, however, a significant normal field through the y = 0 plane appears, and the normal 
field through the z = 0 plane shows signs of becoming asymmetric. Later on (t = 6.0), the By 
pattern at the y = 0 plane becomes stronger, while the pair of Bz patterns at the outer boundaries 
of the two tubes become dominant while the pair at the inner boundaries between the tubes 
gradually merge and cancel. Finally, at t = 10, the By field through they = 0 plane consists of two 
strong, equal and opposite components at the edges of the merged flux tube. The Bz patterns in 
the right panel of the figure generally resemble those of By on the left at t = 10, but Bz still shows 
significantly more structure than By. Nevertheless, the trend is clear: the two flux tubes are 
merging to form a single coherent tube, with the twist concentrated at the outside of the merged 
system. 

We find that magnetic helicity is conserved throughout the evolution; this is demonstrated 
below. Consequently, at least for this simplest possible case, reconnection indeed transports the 
shear to the boundary of the stressed flux system, leaving the central portion relatively smooth, 
as predicted by the helicity condensation model. Furthermore, we have shown that two coronal 
loops with the same sign of helicity reconnect effectively, as conjectured by Antiochos (2013). 

3.2 Counter-Helicity Case 

To examine the interaction of coronal flux tubes with opposite helicity, we simply reverse the 
sense of rotation at the photosphere for one of the flux tubes, as in the right panel of Figure 1, so 
that the two flux tubes are injected with oppositely signed magnetic helicity. Note that the 
geometry and the boundary conditions are kept the exactly the same as in the co-helicity case 
above; the only difference is that the injected helicity has opposite signs in the two flux tubes. 

For comparison with Figure 2, the evolution of the magnetic field is displayed in Figure 5 for 
times t = 3, 4.5, 6, 10. The field lines shown in white and green are traced from points on the 
bottom plane x = 0, along the line z = 0, withy = ±0.005, ±0.0625, ±0.245; the two blue lines are 
traced from the points x = 0, y = 0, and z = ±0. 125. The left colu mn shows top views of the x = 0 
plane, with the plasma velocity magnitude contoured on the plane. The right column shows side 
views of the field lines, with velocity magnitude plotted on the x = 0 (bottom) and x = 0.9 planes. 
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Figure 5. Evolution of eight field lines at selected times t during the counter-helicity simulation. 
Left: X = 0 plane viewed from above; right: x = 0 (bottom) and x = 0.9 planes. Velocity 
magnitude is color-contoured on the planes. Compare with Figure 2 for the co-helicity case. 


Because the magnetic field lines in the two neighboring flux tubes are twisted in opposite 
directions, the tangential components of the magnetic field at the interface of the two tubes (y = 0 
plane) are parallel, rather than anti-parallel as in the co-helicity case. Therefore, we do not 
expect them to reconnect, no matter how much helicity is injected. The results displayed in 
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Figure 5 confirm this expectation: from early on {t = 3) to the final simulation time {t = 10), the 
helicity in each flux tube continuously increases, and the magnetic field lines become more 
twisted. However, no field lines cross the y = 0 plane between the two flux tubes. This is true 
even for the two lines that are rooted very close to the interface between the flux tubes (y = 
±0.005). There is evidence, however, for reconnection between the twisted flux and the external 
untwisted flux. By t = 10, some of the footpoints inside of the twisting regions on the bottom 
boundary now connect to points outside of the twisting regions on the top boundary. Note also 
that the two blue field lines near, but outside, the twisting regions have changed connectivity, 
indicating that they too have undergone reconnection. The key result is that even at the final 
simulation time, t = 10, when the field lines at the maximum angular rotation speed (Qmax) are 
twisted by nearly three full rotations (Qmax^9.75/2;r = 1.876x9. 75/2;i: = 2.91), there is no sign that 
reconnection has linked the two flux tubes. The helicity transport in this case is minimal, at best, 
and there is certainly no evidence of any helicity cancellation even though the total helicity in the 
corona is zero. 

This result can be seen more definitively by examining the evolution of By and as we did 
in Figure 3. Figure 6 shows the contours of those magnetic tangential components on the 
horizontal plane x = 0.5, at the same times as in Figure 5. Unlike Figure 3, there is almost no By 
flux crossing the y = 0 plane at any time during the simulation. For large twist, the By and B^ 
contours become more complex and appear to extend beyond the initial twisting regions due to 
reconnection with external flux, but there are no signs of merging between the two flux tubes. As 
expected, the twisted flux tubes begin to kink and wrap around each other, which causes flux to 
break the y = 0 plane but does not establish any linkages between the tubes. These results 
confirm the conjecture by Antiochos (2013) that opposite-helicity coronal loops do not reconnect 
and, therefore, coronal helicity cannot simply cancel out via reconnection. Furthermore, this 
counter-helicity simulation does not produce magnetic structure compatible with that in the 
observed corona: there is no buildup of twist (shear) at the boundary of the flux system, and the 
flux in the interior becomes increasingly tangled. 

In order to compare the two cases more quantitatively, we plot the evolution of some key 
global quantities for both simulations in Figure 7: (a) total kinetic energy K\ (b) total magnetic 
energy M\ (c) average magnetic flux density Fy through the y = 0 plane; and (d) total magnetic 
helicity H. At early simulation times, the kinetic energy “rings” with small-amplitude peaks at t 
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Figure 6. Left: color contours of By on the x = 0.75 horizontal plane at selected times t during 
the counter-helicity simulation; right: contours of Bz on the same plane. Compare with Figure 3 
for the co-helicity case. 


« 0.75, 1.75, 2.75, and 3.75. The unit period. At = 1, is precisely the time required for an Alfven 
wave to propagate from one photosphere to the other, or from either to the loop apex in the 
corona and back. These peaks evidently result from constructive interference at the loop apexes 
of the waves launched from the photospheres, with a time delay At « 0.25 before the first peak; 
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Figure 7. Global quantities throughout the co-helicity (solid curves) and counter-helicity 
(dashed curves) cases: (a) kinetic energy, K; (b) magnetic energy, M; (c) average magnetic flux 
density through they = 0 plane, Fy; and (d) magnetic helicity, H. 


t = 0.25 is the time at which the imposed photospheric flows attain their maximum acceleration. 
It is evident that over the first three full periods, up to t « 3.75, the kinetic and magnetic energies 
of the two cases are indistinguishable. After this time, the counter-helicity case begins to undergo 
a kink instability, so its kinetic energy starts to rise, whereas its magnetic energy levels off. The 
co-helicity case, on the other hand, continues its early-time behavior until t » 5.5, at which time 
it undergoes a burst of reconnection that sharply increases the kinetic energy and reduces the 
magnetic energy by time t » 6. We note from panel (c) that \By\ increases dramatically over this 
interval as a result of the reconnection. 
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Figure 7(d) shows the time variation of the total magnetic helicity in the two cases. The 
points correspond to values measured at discrete times during the evolution by numerically 
integrating the helicity in the coronal volume, whereas the curves indicate the values expected by 
analytically integrating the helicity injected at the boundary. Of course, for the counter-helicity 
case the total helicity should vanish, and this is exactly what we find. More important, we see 
that the co-helicity simulation accumulates helicity at the peak rate AxdHJdt = 2.4x10”^ after the 
ramp-up phase, so that even though it exhibits copious reconnection, the helicity is conserved 
very accurately. This result is critical: it verifies that our MHD simulations are capturing the 
effects of helicity-conserving reconnection, which is expected to occur in the solar corona, rather 
than dissipating helicity via resistive diffusion of magnetic flux. 

4. INTERACTIONS OF SEVERAL TWISTED FLUX TUBES 

Having established that like-helicity loops within a unipolar region reconnect pairwise, we 
now proceed to test the further conjecture by Antiochos (2013) that the reconnection will 
progress to ever-larger scales, until the twist attains the greatest spatial extent that is available. 
On the Sun, the largest scale is the polarity inversion line (PIT) at the boundary of the unipolar 
region, where the parallel twist fields inside and outside the PIL accumulate and strengthen. Our 
uniform-field configuration has no PIL; rather, we model the effect of the finite extent of a 
unipolar region by limiting the area over which the rotational flows are imposed. In particular, 
we adopt the simplest possible arrangement of a central rotational cell surrounded by six other, 
identical rotations arrayed in a close-packed, regular hexagonal pattern on each photosphere. 
This minimally complex model enables us to examine interactions between multiple like-helicity 
twisted flux tubes within a unipolar region encompassing numerous convection cells. 

The simulation setup and early evolution of this system are shown in Figure 8. We used the 
same vertical extent of the domain as before, but increased the horizontal extent to accommodate 
the larger number of rotations. The resulting domain size is [0,l]x[-0.75,+0.75]x[-0.75,+0.75], 
which, if uniformly gridded, would have been spanned by a 256x384 grid. Paired rotations at 
the bottom (x = 0) and top (x = 1) planes were positioned at identical horizontal centers (yc,Zc), as 
before; for the pattern shown, they were set at yc = 0, Zc = 0; yc = ±2ao, Zc = 0; and yc = ±ao, Zc = 
±V3ao, with ao = 0.125. Vertical planes tangent to each pair of flow circles separate the initially 
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Figure 8. Early evolution (cycle n = 0.5, time t = 1.675) of the simulation with seven rotations in 
a hexagonal pattern. Six field lines with color contours of velocity magnitude on the (a) x = 0 
plane viewed from above and (b) x = 0 (bottom) and x = 0.9 planes. Contours of (c) By and (d) 
Bz on the x = 0.75 horizontal plane. Contours of (e) By on the y = 0 and (f) B^ on the z = 0 
vertical planes. The By and Bz images were enhanced to improve the visibility of the weak fields. 


straight flux tubes that become twisted, with the outer edges of the tubes just touching at those 
planes, as in our earlier experiments with two rotations. 

The rotational flows are smoothly ramped up from zero, to maximum speed, and then back 
down to zero, over twelve cycles {n = 1, 2, 12) each having duration At = 27i/|Q|max = 3.35. 

This procedure induces a maximum twist of one full turn within each flux tube over each cycle. 
The results in Figure 8 were obtained at the midpoint of the first cycle, n = 0.5 (t = 1.675). At 
this twist, the imposed footpoint motions have created seven distinct, separated flux tubes in the 
volume having up to one-half turn of twist each. These features are illustrated by the six field 
lines, one originating in each of the outer flux tubes, drawn in Figures 8(a) and 8(b), and by the 
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views from above of the color-contoured horizontal field components, By and Bz, on the x = 0.75 
plane in Figures 8(c) and 8(d), respectively. On the tangent planes (e.g., y or z) at the boundaries 
between flux tubes, the normal field (By or B^) is zero to a high degree of approximation during 
this early evolution, and negligible flux crosses the planes to link the tubes to one another. 
However, the y = 0 and z = 0 planes are crossed by the growing twist field within individual flux 
tubes, as shown in the color-contoured views of By and Bz, respectively, in Figures 8(e) and 8(f). 
The y = 0 plane shown in (e) corresponds to a horizontal line when viewed from the top in (c,d). 
This line cuts through only the central rotation and contains one flux tube. In contrast, the z = 0 
plane shown in (f) corresponds to a vertical line when viewed from the top in (c,d), and cuts 
through three aligned rotations (including the central one) and contains three flux tubes. At this 
early time, therefore, the magnetic flux densities through the y = 0 and z = 0 planes are 
substantially different. 

In our first experiment described below, we maintain the pattern of seven rotations fixed in 
place as shown in Figure 8, through all eleven successive cycles of ramping up and then ramping 
down the flows. At sufficiently large twist or late times t, we expect the seven flux tubes to 
expand and come into contact, interacting and reconnecting with each other as occurred in our 
previous two-tube, co-helicity simulation. New linkages should form across both the y = 0 and z 
= 0 planes. If the majority of the injected twist flux thereby is transported by the reconnection to 
the outer boundary of the flux region as predicted by Antiochos (2013), furthermore, then the 
entire region should evolve toward a single twisted flux tube, as we observed previously for two 
rotations. This implies that the qualitatively and quantitatively different flux densities shown in 
Figures 8(e,f) should eventually converge and begin to resemble one another. That this indeed 
occurs is demonstrated in the next subsection. 

Of course, the Sun does not exhibit a fixed spatial pattern of surface rotations; individual 
granules and supergranules wax and wane, and the global patterns of their motions are constantly 
shifting. We take into account this stochastic aspect of the solar convection in our second 
experiment described below, to demonstrate that the helicity condensation model is robust under 
this added complexity. At each transition between successive cycles (i.e., at times t = 3.35n for n 
= 1,2, . . ., 1 1), the pattern of seven rotations is randomly shifted and rotated as described below. 
The resultant tangling of neighboring flux tubes redistributes the twist field independently of the 
reconnection that is the sole mode of flux transfer in the case where the rotation pattern is fixed. 
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Nevertheless, we find that the late-time behavior of the large-scale twist field is very similar for 
our fixed and randomized patterns. This result is shown in the second subsection to follow. 

4.1 Fixed Rotation Pattern 

We first consider the case in which the pattern of rotations is spatially fixed, at the location 
whose early evolution was shown in Figure 8. The flow amplitudes are ramped up and down 
cyclically in time, to enable direct comparisons with the corresponding results for the spatially 
randomized pattern of rotations. Figure 9 shows the later evolution of six field lines drawn from 
the same positions as those shown in Figures 8(a,b), at the midpoint of selected cycles {n = 1.5, 
4.5, 7.5, and 11.5) during the simulation. As in the two-tube cases presented earlier, the 
footpoint positions of the field lines are spatially fixed: they are located within the flow annuli of 
the flux tubes, but they do not follow the flow, so in general the field lines are different between 
any two frames. The field-line color uniquely identifies only the starting position from which the 
line is traced. We observe an overall, but non-monotonic, increase in the angle of rotation of the 
field lines about the whole pattern of flows. At n = 1.5, some field lines begin within one flux 
tube at the bottom of the domain (x = 0), but end within a nearby portion of the counter- 
clockwise neighboring flux tube at the top (x = 1), indicating that they already have reconnected 
at least once. The remote end positions of some field lines reaches to the far portion of the 
nearest neighboring flux tube at n = 4.5, and beyond into the second-nearest neighboring tube at 
n = 1 1 .5; but in between, at n = 7.5, the apparent rotation is smaller, more like that at n = 1.5 than 
at the later cycles. Clearly, field lines rooted in the flow annuli endure multiple reco nn ections 
during the evolution. This process transfers magnetic twist flux throughout the volume of axial 
flux linking the surface regions where the rotational flows are imposed, as well as into the 
surrounding line-tied field (note the dark blue field line at n = 11.5). The twist aggregates 
around the perimeter of the flow pattern: no field line ends within the central flux tube at any 
time shown, and only the green field line at n = 1 1.5 threads tangentially through its periphery. 

Figure 10 shows views from above of the contours of By (left panels) and (right panels) on 
the horizontal plane x = 0.75, at the same times as the field lines drawn in Figure 9. The images 
reveal the gradual distortion, partial cancellation, and large-scale merger of the individual flux 
tubes that were so clearly evident (after enhancement to emphasize weak fields) in Figure 8(c,d). 
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Figure 9. Evolution of six field lines at selected cycles n (times t = 3.3 5n) during the fixed- 
pattern simulation. Left: x = 0 plane viewed from above; right: x = 0 (bottom) and x = 0.9 
planes. Velocity magnitude is color-contoured on the planes. 


At the final cycle shown {n = 1 1.5), the twist components of the field have “condensed” into two 
strong bands of opposite sign, resembling a single flux tube encompassing the entire region of 
flow. This evolution closely matches that found in our co-helicity simulation (Figure 3), but for 
several more rotational flows (seven vs. two) and, therefore, over a larger range of flux scales. 
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Figure 10. Left: contours of By on the x = 0.75 horizontal plane at selected cycles n (times t = 
3.35n) during the fixed-pattern simulation. Right: contours ofB^ on the same plane. 


As occurred in that previous simulation, the pattern of B^ contours (right) at the last cycle {n = 
11.5) in Figure 10, if rotated clockwise by 90°, would resemble closely that of By (left) at the 
same time in the figure. 
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Figure 11. Left: contours of By on they = 0 vertical plane at selected cycles n (times t = 3.35n) 
during the fixed-pattern simulation. Right: contours of Bz on the z = 0 vertical plane. 

Similar processes of distortion, cancellation, merging, and strengthening of the field patterns 
are observed in Figure 11, which shows the normal By and Bz field eomponents through thej = 0 
and z = 0 planes on the left and right, respectively. The first signs of formation of outer bands of 
By flux, due to linkages across the y = 0 plane established by recormeetions between flux tubes, 
already are evident at n = 1.5. These linkages compete in strength with the twist field indueed 
direetly within the eentral flux tube, whereas only the latter is evident at n = 0.5 in Figure 8(e). 
Meanwhile, the Bz fluxes across the z = 0 plane are strengthening at heights x =; 0.5 at n = 1.5 in 
Figure 11, where the flux tubes beeome erowded together as they expand laterally in response to 
the increased magnetie pressure. As time progresses, the irmer patterns of By and Bz twist field 
associated with individual flux tubes weaken and fade, giving way to a single pair of opposite- 
polarity bands of flux at the perimeter of the region of rotational flow. This late-time structure 
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again resembles a single flux tube encompassing the entire flux system (cf. Figure 4 for the co- 
helicity simulation). As anticipated, the initially very different y and z flux densities have 
changed to resemble one another quite closely at late times in the simulation. 

It is demonstrated below that the magnetic helicity is conserved throughout this evolution. 
Thus, for the case of several rotational flows in a spatially fixed pattern, we again find that 
reconnection transports the magnetic shear to the boundary of the flux system while leaving the 
central portion relatively smooth. This result is in accord with the predictions of the helicity 
condensation model (Antiochos 2013). 

4.2 Randomized Rotation Pattern 

In order to investigate the effect of a constantly shifting pattern of flows, as occurs in the 
Sun’s supergranulation, we repeated the simulation above while adding a random lateral 
translation and angular rotation of the entire hexagonal pattern at the beginning of each flow 
cycle after the first. The original pattern obviously is invariant to a rotation of 60° about its 
center point. In addition, a lateral translation of one flow diameter, do = 2ao, along any of its 
three symmetry lines separated by 60°, places four of the shifted rotations on top of four of the 
original rotations; thus, if the pattern were unbounded rather than bounded in extent, it would be 
invariant to such a displacement. Therefore, beginning with the hexagon centered at (yh,Zh) = 
(0,0) and having a rotation angle = 0 relative to the y axis of the coordinate system, as shown 
in Figure 8, we generated independent sequences of random numbers to set new center 
coordinates y\!,z\! e [-ao,+ao] and a rotation angle (j\! e [-30°, +30°] for each of the ensuing 
eleven flow cycles. Table 1 lists the resulting set of parameters, which were used to impose the 
same shifted and rotated pattern at both the bottom (x = 0) and top (x = 1) photospheres in each 
cycle. This randomization process broadly distributed the twist imparted to the magnetic field 
across the footprint of the entire hexagonal pattern of motions, rather than concentrating it in 
seven discrete, fixed flux tubes. However, as will be seen below, the late- time changes to the 
magnetic field of the randomized-pattem simulation relative to the fixed-pattern case are not 
substantial, having a more quantitative than qualitative character. 
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n 

Th 

Zh 


1 

0.0 

0.0 

0° 

2 

+0.1020 

+0.0745 

-11° 

3 

-0.0130 

-0.0155 

+10° 

4 

+0.0320 

-0.0480 

+22° 

5 

+0.1160 

-0.0405 

+23° 

6 

+0.0920 

-0.0990 

-6° 

7 

+0.1070 

-0.1100 

+27° 

8 

-0.0685 

+0.0040 

O 

O 

(N 

1 

9 

+0.0065 

+0.0115 

-15° 

10 

+0.0760 

-0.0110 

O 

(N 

1 

11 

+0.1035 

+0.0645 

-10° 

12 

+0.1200 

-0.0865 

+19° 


Table 1: Parameters yh, Zh, and ^ of the randomized pattern of rotations vs. convection cycle n. 

Six selected magnetic field lines from the randomized-pattem simulation are shown in Figure 
12 at the same selected cycles {n = 1.5, 4.5, 7.5, and 11.5) as were shown for the fixed-pattern 
simulation in Figure 9. In this case, the footpoint locations are not fixed in space for all time, as 
they were before, but instead are translated and rotated to follow the randomization of the flow 
pattern. The resultant shifting of the pattern during this sequence is easily seen in the contour 
shading of velocity magnitude on the x = 0 and x = 0.9 planes. Less obvious, but also apparent, 
is the effect of the pattern shifting on the twist associated with individual field lines: the angular 
rotation of the lines in the early cycles {n= 1.5 and 4.5) is, on average, smaller for the random 
pattern than it is for the fixed pattern at the same time. This is unsurprising, given the loss of 
coherence of the twist applied to particular magnetic flux surfaces when the pattern of flows is 
periodically shifted, compared to when it is fixed. Nevertheless, the field line rotation in the late 
cycles {n = 7.5 and 1 1.5) is quite large and rather similar for the random and fixed patterns, and 
the twist aggregates near the perimeter of the random flow pattern as it did for the fixed. 

These trends are seen more clearly in views of the By and contours from above (Figure 13) 
and from the sides (Figure 14). The coherent small-scale structure of individual flux tubes in the 
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Figure 12. Evolution of six field lines at selected cycles n (times t = 3.3 5n) during the 
randomized-pattern simulation. Left: x = 0 plane viewed from above; right: x = 0 (bottom) and x 
= 0.9 planes. Velocity magnitude is color-contoured on the planes. Compare with Figure 9 for 
the fixed-pattern simulation. 

fixed-pattern simulation (Figures 10 and 11) fragments and fades noticeably more rapidly in the 
randomized-pattern simulation (Figures 13 and 14) throughout the early and middle cycles {n = 
1.5, 4.5, and 7.5). Consequently, the emerging large-scale structure, which resembles a single. 
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Figure 13. Left: contours of By on the x = 0.75 horizontal plane at selected cycles n (times t = 
3.35n) during the randomized-pattern simulation. Right: contours of Bz on the same plane. 
Compare with Figure 10 for the fixed-pattern simulation. 

twisted flux tube encompassing the entire region of rotational flow, is more obviously apparent 
during the middle cycles {n = 4.5 and 7.5) of the random pattern than of the fixed. By the end of 
the simulations {n = 11.5), however, the large-scale structure has become dominant and the two 
late-time configurations are essentially indistinguishable. 
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Figure 14. Left: contours of By on they = 0 vertical plane at selected cycles n (times t = 3.35n) 
during the randomized-pattern simulation. Right: contours of Bz on the z = 0 vertical plane. 
Compare with Figure 11 for the fixed-pattern simulation. 

A more quantitative comparison of the two cases is shown in plots of the evolution of key 
global quantities. Figure 15 displays histories of (a) total kinetic energy K, (b) total magnetic 
energy M, and (c) total magnetic helicity H. The kinetic and magnetic energies for the fixed and 
randomized patterns are identical, of course, over the first cycle (to time t = 3.35), and thereafter 
exhibit qualitatively similar behaviors. The kinetic energies in Figure 15(a) attain maxima at 
approximately the midpoints of the cycles (at times t « 3.35n+1.675, n > 1), followed soon 
thereafter by peaks in the magnetic energies, and minima at approximately the endpoints of the 
cycles (at times t « 3.35n, n > 2). Decreases in the magnetic energies in Figure 15(b) indicate 
the occurrence of reconnections between flux tubes in the case with fixed rotations (solid curve)', 
both reconnections between, and rotation-induced untwisting of, flux tubes occur in the case with 
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Figure 15. Global quantities throughout the fixed-pattern (solid curves; circles) and 
randomized-pattern (dashed curves; squares) simulations: (a) kinetic energy, K; (b) magnetic 
energy, M; and (c) magnetic helicity, H. 


randomized rotations (dashed curve). When these processes come into play, early in the second 
cycle, the overall trend lines of increasing magnetic energy adopt gentler slopes and the maxima 
in the kinetic energy increase substantially. 

Figure 15(c) shows the evolution of the total magnetic helicity in the two cases. The points 
correspond to values measured at discrete times during the evolution by numerically integrating 
the helicity in the coronal volume, whereas the curve indicates the co mm on value expected by 
analytically integrating the helicity injected at the boundary. The two sets of data points overlie 
each other, because randomizing the flow pattern does not alter the pace of helicity injection, 
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Figure 16. Global quantities throughout the (a) fixed-pattern and (b) randomized-pattern 
simulations: average magnetic flux densities through the y = 0 plane, Fy (solid curves), and the z 
= 0 plane, F^ (dashed curves). 

whose average rate is O.SxXAxdHJdt = 4.2x10 . As in the two-tube co-helieity case discussed 
previously, even though both the fixed- and randomized-pattem simulations exhibit copious 
reconnection, the helicity is conserved very accurately throughout. 

Figure 16 displays the average magnetic flux densities, Fy through the y = 0 plane and F^ 
through the z = 0 plane, for the (a) fixed-pattern and (b) randomized-pattem simulations. At very 
early times, these quantities reflect the formation of the individual twisted flux tubes in response 
to the rotational motions imposed at the photospheres. As noted in conjunction with Figure 8, 
the initial orientation of the hexagonal flow pattern is such that only one flux tube cuts through 
the y = 0 plane while three tubes cut through the z = 0 plane. This imparts the difference of a 
factor of three in the initial slopes of the Fy (solid) and F^ (dashed) curves in Figure 16. During 
the second cycle of flows - at its midpoint for the fixed pattern in 16(a) and at its beginning for 
the randomized pattern in 16(b) - the trend in F^ moderates strongly as reconnections between 
flux tubes cancel flux crossing through the z = 0 plane and, in the case of the randomized pattern, 
the initial alignment of the three rotational flows along the z = 0 line is lost. In contrast, the 
initial trend in Fy is maintained, with small-amplitude transient oscillations superimposed on the 
trend line. At late times in the fixed-pattern simulation in 16(a), the two curves converge, 
signaling the final loss of memory of the initially anisotropic configuration during the approach 
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to an isotropic state resembling a single twisted flux tube. This convergence and loss of memory 
occur quite early in the randomized-pattem simulation in 16(b), due to the periodically shifting 
pattern of flows. Thereafter, the two curves trend together and oscillate about each other. At the 
end of the simulations, the flux densities for the fixed and randomized patterns are essentially 
identical along both the y and z directions and to each other. 

These results demonstrate that the long-term averaged behavior of coronal flux systems 
subject to helicity injection, transport by reconnection, and condensation at the boundary is not 
sensitive to the steady or unsteady character of the pattern of surface rotations. Regardless of 
whether the pattern of flows is fixed or varies randomly, over time the twist field concentrates at 
the perimeter of the region of imposed flows. Furthermore, at identical rates of helicity injection, 
the accumulated twist flux becomes independent of the fixed or randomized nature of the flows 
once the memory of the particular initial and boundary conditions imposed has been lost. This is 
an important additional verification of the robustness of the helicity condensation model. 

5. DISCUSSION 

The results presented in this paper yield several important conclusions on the structure and 
dynamics of the coronal magnetic field. Clearly, these first, exploratory simulations do not 
include all of the observed properties of the coronal field, but they do capture some essential 
features. Perhaps the most essential feature is that the simulations conserve magnetic helicity to 
an excellent approximation, even though there is clearly a great deal of reconnection occurring in 
the coronal domain. We note that in Figure 15c, the points corresponding to the numerically 
measured helicity fall almost exactly on the curve corresponding to the analytic expected value 
(within « 1%). Thus, the evolution determined by our ARMS simulations is dominated by true 
reconnection, rather than diffusion. There is good reason to expect, therefore, that our results are 
robust in that they are unlikely to change when scaled up to a Lundquist number appropriate to 
the real corona, which is orders of magnitude larger than in any possible simulation. One 
uncertainty is whether reconnection occurs as readily in the corona as in our simulations. A 
number of studies, however, have shown that complex photospheric flows with stagnation points 
produce current structures in the corona with widths that decrease exponentially over time (e.g., 
van Ballegooijen 1986; Antiochos & Dahlburg 1997). Consequently, even for coronal Lundquist 
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numbers, we expect that the current structures will thin down to the dissipation scale and start to 
reconnect in only a few rotations. 

A key conclusion of our simulations is that the helicity evolution of the corona is highly 
constrained. We find that the corona clearly does not evolve toward the Taylor (1986) state of a 
linear force-free field. For the simulations with two flux tubes, a Taylor evolution implies the 
complete cancellation of all magnetic twist in the counter-helicity case and the spreading of twist 
uniformly throughout the domain in the co-helicity case. In the counter-helicity case, we find 
instead that the magnetic twist remains concentrated in the two separate flux tubes, with little 
evidence for reconnection even after the tubes kink substantially. In the co-helicity case, the twist 
evolves via reconnection to merge the two flux tubes, but it does not propagate significantly into 
the surrounding field. Note that even the simulations with multiple flux tubes, which involve a 
great deal of reconnection, do not lead to a minimum-energy Taylor state. 

These outcomes occur because line-tying at the photosphere introduces strong dynamical 
constraints on the coronal evolution. As discussed in Antiochos et al. (2002), reconnection 
requires current sheets, but these cannot form freely throughout the corona. They are likely to 
form only where the field has topological singularities, such as null points and separatrices, or 
where the driving motions have topological singularities, such as stagnation points. For a bipolar 
coronal region with a topologically smooth potential field that is continuously driven by 
photospheric motions, we find that the helicity condenses at the largest system scale, in 
agreement with turbulence theory, rather than evolving toward a Taylor state. 

Even in our simulation with random photospheric motions, which inherently contains many 
stagnation points and consequently numerous current sheets, the system evolves via helicity 
condensation rather than to a linear force-free state. Due to the low beta of our simulation the 
final state is, indeed, force-free to a good approximation, but it is far from a linear force-free 
state. Also, for the random driving we expect that, in addition to simply helicity, higher-order 
topological features such as braiding are injected into the coronal field. It appears, however, that 
all of these higher-order topologies are destroyed by reconnection and only helicity remains, in 
agreement with other calculations (e.g., Pontin et al. 2011). 

The conclusion that coronal magnetic helicity evolves via condensation to the largest scale is 
the most important result of our simulations. It has far-reaching implications for solar 
observations; in particular, it explains why coronal loops appear laminar while shear collects at 
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polarity inversion lines. The calculations presented here provide strong initial support for the 
conjectures developed in Antiochos (2013), but much more observational and theoretical work 
remains to be done in order to confirm or refute the helicity condensation model. 
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